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Abstract. Based on ideas by Woodward et al., a subgrid scale model that 
is applicable to highly compressible turbulence is presented. Applying the sub- 
grid scale model in large eddy simulations of forced supersonic turbulence, the 
bottleneck effect is largely reduced and, thereby, approximate scaling laws can 
be obtained at relatively low numerical resolution. In agreement with previous 
results from PPM simulations without explicit subgrid scale model, it is found 
that the energy spectrum function for the velocity field with fractional density- 
weighing, p 1 / 3 !*, varies substantially with the forcing, at least for the decade of 
wavenumbers next to the energy-containing range. Consequently, if universal 
scaling of compressible turbulence exists, it can be found on length scales much 
smaller than the forcing scale only. 



1. Introduction 

Large eddy simulations (LES) are of great utility to engineers and atmospheric 
scientists, but not commonly used in computational astrophysics. While terres- 
trial turbulence is incompressible or weakly compressible in most cases, astro- 
physicists often deal with highly compressible turbulence. Since the stability of 
numerical solvers for the equations of compressible gas dynamics requires energy 
dissipation that is intrinsic to the numerical scheme, gas-dynamical simulations 
at high Reynolds numbers are considered as implicit large eddy simulations. 
However, it was s hown that numerical schemes s uch as the piecewise parabolic 
method (PPM) of IColella and Woodwardl ( 19841 ) entail undesired properties of 



the numerical solutions such as the bottleneck effect, i. e., an unphysical en- 
hanc ement of s pectral power in the high -wave number range (e. g.. lKritsuk et al 



20071 ). This led I Woodward et all (|2006l ) to the idea to couple an explicit subgrid- 



scale (SGS) model to PPM, and they applied the method to decaying transonic 
turbulence. 

In this article, I briefly outline an SGS model that is valid in the supersonic 
regime. As a first application, I present results from LES of forced supersonic 
isothermal turbulence with root mean square Mach number between 5 and 6, 
which can be interprete d as an idealized model for cold turbule nt gas in the 
interstellar medium (see lKritsuk et ai1l2007l ; iFederrath et aZ.| [2009) . Computing 



turbulence energy spectra from these LES, I find that the SGS model largely 
reduces the bottleneck effect in comparison to plain PPM simulations. The 
results shed light on the question of the universality of compressible turbulence, 
which is important for the theory of star formation in molecular clouds (see, for 
example, lElmegreen and Scald 12003 ) . 
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2. A Subgrid Scale Model for Supersonic Turbulence 



ISchmidt et al\ ( 20061 ) formulated an SGS model based on the energy contained in 
numerically unresolved velocity fluctuations, K sgs , to compute turbulent mixing 
properties. For example, this model improved the treatment of turbulent com- 
bustion in thermonuclear supernova simulations. The major modeling challenge 
is to compute a stress tensor, r sgs , that accounts for the interaction between nu- 
merically resolved scales and the subgrid scales below the grid resolution. The 
divergence of this stress tensor, V -Tsgs, enters as an additional force term in the 
Euler equation for the conservation of momentum. While the diagonal part of 
T sgs introduces quasi- viscous effects analogous to the viscous term in the Navier- 
Stokes equation, the diagonal part corresponds to an SGS turbulence pressure 
that alters the effective equ ation of state in the highly compressible case. 

Woodward et al\ (|2006l ) proposed a non-linear closure for r sgs , which is sec- 



ond order in the trace-free part of the velocity derivative V ® u rather than first 
order in the trace-free rate of strain (as in the standard eddy- viscosity closure). 
I invest igated a closure that is a combination of the linear eddy- viscosity closure 
used in lSchmidt et al. ( 20061 ) and a non- linear term that is constructed from the 



complete velocity derivative: 



Tij = CxAp^KlgStj - 2C 2 K sgs ^ 1 ^ - |(1 - C 2 )if sgs <%, (1) 



where A is the size of the grid cells, p is the mass density, | V<g>u| := (2-Uj^Uj^) 1 / 2 , 
and S*j is trace- free symmetric part of Uij. The SGS turbulence energy density 

-Ksgs := —\ T a is given by a dynamical equ ation that compleme nts the Euler 
equations for compressible gas dynamics (see lSchmidt et al.ll2006l ). 



In orde r to verify the applicab i lity o f this closure to isotropic supersonic 
turbulence, ISchmidt and Federrathl ( 2009h filtered numerical data from 1024 3 

up 



PPM simulations (jFederrath et aiJ l2009) on a length scale within the approx- 
imate inertial subrange. Comparing the explicit computation of the rate of 
energy transfer, £ = TijUij, from length scales greater than the filter length to 
smaller length scales to the predictions following from various closures for r«, 
we found that the closure ([T]) yields the best correlation. The closure validation 
as well as further details of the SGS model will be described in a separate article 



(ISchmidt and Federrathl 120091 ) . 

As an illustration, Fig.[T]shows a visualization of K sgs prepared from an LES 
with 512 3 grid cells. In this simula tion, purely solenoid al forcing was applied (the 
forcing algorithm is described in Schmidt et al. 2009). In the reddish regions, 



K sgs is higher than the average value, while it is lower in bluish regions. It turns 
out that K sgs is fairly well correlated with the local denstro phy ^ | V x fpy 2 tt) |, 
an indicator of compressible turbulent velocity fluctuations ( Kritsuk et al. 20071 ). 



This reflects the locality of turbulence energy transport from resolved to unre- 
solved scales, although the distribution of K sgs is smoother in comparison to the 
local denstrophy, because of the diffusive properties of the SGS dynamics. 
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Figure 1. Visualization of the SGS turbulence energy density K sgs in a 512 3 
LES with purely solenoidal forcing. 



3. Turbulence Energy Spectra 



It is not obvious which quantity should be considered to specify spectral prop- 
erties of highly compressible turbulence. Numerical simulations of supersonic 
turbulence driven by a stirring force showed that the spectrum functions E{k) of 
kinetic energy per unit mass are significant ly stiffer than a Kolmogorov spectrum 
with slope —5/3 (e. g. jKritsuk et al. 20071). and indications for a dep endence on 
the forcing were found (jSchmidt et a{j|200<l iFederrath efaR 120091 ). However, 
because of the huge fluctuations of the mass density, it seems natural to apply 
some kind of density weighing. A generalized turbulence energy spectrum func- 
tion with fractional density weighing, E q (k), is given by integrals over spherical 
surfaces of radius k in Fourier space, 



EJk) 



\( P i u ) 2 (k)\k 2 dn h , 



(2) 



where p q u is the Fourier transform of p q u. Since the kinetic energy density 
is given by \pu 2 , calculating spectrum functions Ei/ 2 {k) is an option that sug- 
gests itself. Numerical results indicate that the slope of these spectrum functions 
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tend to be natter than -5/3 (|Kritsuk et a/.ll2007l : Ifbderrath et al.ll2009l ). On the 
other hand, iLighthilll (| 19551 ) recognized the significance of the rate of dissipation 
per unit volume for compressible turbulence, which is proportional to the mass 
density times the third power of the velocity fluctuation on a given length scale. 
A statistically constant dissipation rate in the inertial subrange implies that 
the spectrum function Ei^(k) for the density-weighted velocity field p 1 / 3 u has 
the slope —5/3. Following this reasoning, iKritsuk et all ( 20071 ) demonstrated 
that the spectrum function Ei/^{h) computed from a simulation of supersonic 
isothermal turbulence, indeed, is close to a Kolmogorov spectrum in the iner- 
tial subrange. In the following, I present some results concerning the spectral 
properties of supersonic turbulence computed by means of LES. 

Let us first comp are spectra from LES v s. plain PPM simulations on grids 
with 256 3 cells. As in lFederrath et all (2X109), I consider the two limiting cases 
of purely solenoidal (divergence-free) and compressive (rotation-free) stochas- 
tic force fields. Time-averaged compensated spectrum functions, k 5 ^ 3 Ei/ 3 (k), 
calculated from these simulations are plotted in Fig. [2j In contrast to the LES 
spectra (right plot), bumps caused by the bottleneck effect at high wavenumbers 
noticeably flatten the spectra obtained from the corresponding PPM simulations 
without SGS model (left plot). Since the range of wave numbers, for which the 
spectrum function drops rapidly due to viscous effects, is about the same in both 
cases, the action of the SGS turbulence stress does not significantly degrade the 
dynamical range of PPM, which is an important feature. 




Figure 2. Compensated spectrum functions, k 5 ^ 3 Ei/ 3 (k), for the mass- 
weighted velocity field p x ^u in plain PPM simulations (left) and LES (right) 
on 256 3 grid. The spectra for solenoidal and compressive forcing are shown in 
blue and red, respectively. The dot-dashed lines show transversal components, 
and the dashed lines show longitudinal components of Ei/ 3 (k). 

Fig. [2] suggests that in the intermediate wavenumber range Ei/${k) is close 
to a Kolmogorov spectrum for solenoidal forcing, while it is much stiffer in the 
case of compressive forcing. One can also see that the contributions from the 
transversal and the longitudinal components of p 1 / 3 u to the spectrum functions 
are about the same regardless of the forcing. Spectra for different mixtures of 
solenoidal and compressive forcing modes are shown in the left plot of Fig. [3J 
The slopes are close to —5/3 if the forcing is mainly solenoidal, but the spectra 
become substantially steeper if the forcing is dominated by compressive modes. 
To determine the slopes of the spectrum functions in between the peaks at low 
wave numbers and the viscous cutoff at high wavenumbers, I performed LES also 
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on 512 3 grids. Power-law fits E-^j^k) oc k P 1 / 3 yield the spectral indices /3i/ 3 = 
1.61 and 2.16 for solenoidal and compressive forcing, respectively (see the right 
plot in Fig. [3]). In the latter case, there is a pronounced bend at 2irk/L ~ 10 (L 
is the forcing length scale), and the spectrum function almost flattens to a slope 
of —5/3 at higher waver numbers. I discuss this feature, which does not emerge 
until a resolution of at least 512 3 is used, in the concluding Section. In the case 
of solenoidal forcing, the spectral index slightly decreases at higher resolution, 
while the shape of the spectrum functions does not change. Nevertheless, since 
density-weight ed velocity scalings tend to be very sensitive on the numerical 
resolution (see lSchmidt et a71l2008l ). it is encouraging that the spectral indices 
inferred from the LES agree quite well with the val ues Pi/3 = 1^64 (solenoidal 
forcing) and 2.10 (compressive forcing) reported by IFederrath et al\ (|2009l ) for 
the 1024 3 plain PPM data. 
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Figure 3. Left plot: Results from 256 3 LES, whe re the weighing parameter 
C of the forcing (see eqn. (8) in lSchmidt et a/.ll2009t ) ranges from ( = I (purely 
solenoidal, shown in blue) over 2/3, 1/3, and 1/6 to (purely compressive, 
shown in red). Right plot: Compensated spectrum functions for 512 3 LES 
with solenoidal (blue) and compressive (right) forcing. 



4. Discussion and Conclusion 

The bottleneck effect, which distorts turbulence energy spectra toward high 
wavenumbers in PPM simulations without explicit SGS model, is mostly com- 
pensated by the SGS turbulence stress in LES. This property is clearly favor- 
able for the determination of scaling laws from the spectra. Moreover, the SGS 
model for compressible turbulence presented in this article will be valuable in 
simulations of complex astrophysical flows, where a large range of scales covers 
energy-injecting physical processes, while only a small portion of the inertial 
subrange can be r esolved. A first exa mple for such an applications has recently 
been presented bv lMaier et al\ ( 20091 ) 



Based on new data for hydr odynamical as well as magnetohydrodynamical 
turbulence. iKritsuk et al\ ( 20091 ) put forward the hypothesis that the universal- 



ity of compressible turbulence becomes manifest in the scalings of p^/^u. The 
spectral indices inferred from the 512 3 LES presented in this article, however, 
do not support this hypothesis. There are two interpretations of these results. 
One interpretation presumes that the scalings in the range of wavenumbers next 
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to the energy-containing range reflect genuine inertial-range properties of su- 
personic turbulence. Then the turbulence energy sp ectra clearly exhibit a de- 
pendence of these properties on th e forcing (also see ISchmidt et al\ 120081 ) . As 
proposed by iFederrath et all ( 20091 ). the prominent flattening of the spectrum 



function Ei/ 3 (k) at high wavenumbers that is observed for compressive forc- 
ing (see right plot in Fig. ED might be caused by a transition to the weakly 
compressible regime at the sonic wavenumber, for which the turbulent velocity 
fluctuations are comparable to the speed of sound. On the other hand, it is pos- 
sible that compressible turbulence with universal scaling can only develop on 
length scales much smaller than forcing scale if the forcing is mostly compres- 
sive. In order to settle the question which interpretation is correct, the combined 
influence of varying the Mach number of the flow and choosing different mixing 
ratios of solenoidal and compressive forcing components has to be investigated. 
In this regard, the SGS model presented in this article can help to explore the 
parameter space. 

Acknowledgments. I thank Christoph Federrath for performing the 1024 3 
PPM simulations that served as test cases for the LES described in this article 
and Jens Niemeyer for valuable comments. The SGS model was implemented 
into the Enzo 1.5 code developed by the Laboratory for Computational Astro- 
physics at the University of California in San Diego (htt p://lca.ucsd.edu[ ) 



References 

Colella, P. and Woodward, P. R. (1984). The Piecewise Parabolic Method (PPM) for 
Gas-Dynamical Simulations. Journal of Computational Physics, 54, 174-201. 

Elmcgrccn, B. G. and Scalo, J. (2004). Interstellar Turbulence I: Observations and 
Processes. ARA&A, 42, 211-273. 

Federrath, C, Duval, J., Klessen, R., Schmidt, W., and Mac Low, M. M. (2009). Com- 
paring the statistics of interstellar turbulence in simulations and observations. 
E-print arXiv:0905.1060. 

Kritsuk, A. G, Norman, M. L., Padoan, P., and Wagner, R. (2007). The Statistics of 
Supersonic Isothermal Turbulence. ApJ, 665, 416-431. 

Kritsuk, A. G, Ustyugov, S. D., Norman, M. L., and Padoan, P. (2009). Simulations 
of Supersonic Turbulence in Molecular Clouds: Evidence for a New Universality. 
In N. V. Pogorelov, E. Audit, P. Colella, & G. P. Zank, editor, Astronomical 
Society of the Pacific Conference Series, volume 406 of Astronomical Society of 
the Pacific Conference Series, pages 15 — h 

Lighthill, M. J. (1955). The Effect of Compressibility on Turbulence. In Gas Dynamics 
of Cosmic Clouds, volume 2 of IAU Symposium, pages 121 — h 

Maier, A., Iapichino, L., Schmidt, W., and Niemeyer, J. C. (2009). Adaptively refined 
large eddy simulations of clusters. E-print arXiv:0909.1800. 

Schmidt, W. and Federrath, C. (2009). In preparation. 

Schmidt, W., Niemeyer, J. C, and Hillebrandt (2006). A localised subgrid scale model 
for fluid dynamical simulations in astrophysics. I. Theory and numerical tests. 
Astron. & Astrophys., 450, 265-281. 

Schmidt, W., Federrath, C, and Klessen, R. (2008). Is the Scaling of Supersonic Tur- 
bulence Universal? Phys. Rev. Lett., 101, 194505. 

Schmidt, W., Federrath, C, Hupp, M., Kern, S., and Niemeyer, J. C. (2009). Numerical 
simulations of compressively driven interstellar turbulence. I. Isothermal gas. 
A&A, 494, 127-145. 



LES of Supersonic Turbulence 



7 



Woodward, P. R., Porter, D. H., Anderson, S., Fuchs, T., and Herwig, F. (2006). 
Large-scale simulations of turbulent stellar convection flows and the outlook for 
petascale computation. Journal of Physics Conference Series, 46, 370-384. 



